#include <math.h>

#include <iostream>

#include "Interpolation2D.h"
#include "test.h"

typedef double (*Func)(double, double);

template <class InternalType = double, class ExternalType = float>
void test_subcase(Interpolation2DBase<InternalType, ExternalType> &b,
                  ExternalType xa[], ExternalType ya[], ExternalType **za,
                  int m, int n, Func f) {
  b.construct(xa, ya, za);
  for (int i = 0; i < m; i++) {
    for (int j = 0; j < n; j++) {
      ExternalType x = sqrt(i);
      ExternalType y = sqrt(j);
      // 1e-6 because float has 6 exact numbers
      ExternalType error = (f(x, y) == 0)
                               ? fabs(b.evaluate(x, y) - f(x, y))
                               : fabs(b.evaluate(x, y) - f(x, y)) / f(x, y);
      CHECK(error < 1e-6);
      // printf("x = %.1f, y = %.1f, f = %.2f\n", x, y, b.evaluate(x, y));
    }
  }
}

template <class InternalType = double, class ExternalType = float>
void test_case(int m, int n, Func f, int deg) {
  ExternalType *xa = new ExternalType[m];
  ExternalType *ya = new ExternalType[n];
  ExternalType **za = new ExternalType *[m];
  for (int i = 0; i < m; i++) xa[i] = sqrt(i);
  for (int j = 0; j < n; j++) ya[j] = sqrt(j);
  for (int i = 0; i < m; i++) {
    za[i] = new ExternalType[n];
    for (int j = 0; j < n; j++) za[i][j] = f(xa[i], ya[j]);
  }
  // SUBCASE("Polynomial Interpolation") {
  //   // std::cout << "Polynomial Interpolation" << std::endl;
  //   switch (deg) {
  //     case 1: {
  //       PolInt2D<InternalType, ExternalType, 1> p1 =
  //           PolInt2D<InternalType, ExternalType, 1>(m, n);
  //       test_subcase(p1, xa, ya, za, m, n, f);
  //       break;
  //     }
  //     case 2: {
  //       PolInt2D<InternalType, ExternalType, 2> p2 =
  //           PolInt2D<InternalType, ExternalType, 2>(m, n);
  //       test_subcase(p2, xa, ya, za, m, n, f);
  //       break;
  //     }
  //     default: {
  //       PolInt2D<InternalType, ExternalType, 3> p3 =
  //           PolInt2D<InternalType, ExternalType, 3>(m, n);
  //       test_subcase(p3, xa, ya, za, m, n, f);
  //       break;
  //     }
  //   }
  // }
  SUBCASE("Efficient Linear Interpolation") {
    // std::cout << "Efficient Linear Interpolation" << std::endl;
    LinInt2D<InternalType, ExternalType> p =
        LinInt2D<InternalType, ExternalType>(m, n);
    test_subcase(p, xa, ya, za, m, n, f);
  }
  // SUBCASE("Hermite Interpolation") {
  //   // std::cout << "Hermite Interpolation" << std::endl;
  //   HerInt2D<InternalType, ExternalType> h =
  //       HerInt2D<InternalType, ExternalType>(m, n);
  //   test_subcase(h, xa, ya, za, m, n, f);
  // }
  SUBCASE("Quadratic Interpolation") {
    // std::cout << "Quadratic Interpolation" << std::endl;
    QuaInt2D<InternalType, ExternalType> h =
        QuaInt2D<InternalType, ExternalType>(m, n);
    test_subcase(h, xa, ya, za, m, n, f);
  }
  SUBCASE("Spline Interpolation") {
    // std::cout << "Spline Interpolation" << std::endl;
    SplInt2D<InternalType, ExternalType> s =
        SplInt2D<InternalType, ExternalType>(m, n);
    test_subcase(s, xa, ya, za, m, n, f);
  }
  for (int i = 0; i < m; i++) delete[] za[i];
  delete[] xa, ya, za;
}

TEST_CASE("Test: Bi-Linear Polynomial") {
  int m = 10, n = 10;
  test_case<double, float>(m, n, plinear, 1);
}

TEST_CASE("Test: Bi-Quadratic Polynomial") {
  int m = 10, n = 10;
  test_case<double, float>(m, n, pquadratic, 2);
}

TEST_CASE("Test: Bi-Cubic Polynomial") {
  int m = 10, n = 10;
  test_case<double, float>(m, n, pcubic, 3);
}

TEST_CASE("Test: Trigonometric Function") {
  int m = 10, n = 10;
  test_case<double, float>(m, n, fsin, 3);
}
